Introducción

La encuesta domiciliaria de origen-destino en hogares ATU 2023 se compone de 3 bases de datos:

El Area de Lima Metropolitana esta divido en 1,209 zonas, las cuales no son comparables con las ZAT de la encuesta que hizo la JICA en el 2012. Sin embargo, existe un problema de muestreo ya que solo 641 de estas 1,209 zonas (el 53%) tienen por lo menos un hogar encuestado. Entonces estas zonas no son adaptadas para calcular indicadores a la escala infra-distrital. Como solución, hemos explorado sustituirlas por las ZAT de la encuesta del 2012. Resulta que 344 de las 409 ZAT (84%) tienen por lo menos un hogar encuestado, lo cual es mucho más adaptado. Sin embargo, ante el efecto gráfico de la ausencia de datos en el 16% de las zonas, adoptamos una representación a nivel de distrito.

Nota sobre los factores de expansión

Se manejan varios factores de expansión:

Sin embargo, encontramos una limitación con estos factores de expansión que al parecer no representan bien la distribución de la población dentro de los distritos. Asimismo, lineas abajo comparamos en una tabla la población de los distritos segun proyecciones del INEI para el año 2024 a partir del censo oficial del 2017 y segun la suma de los factores de expansión fexp_per. Como solución, volvemos a calcular nuevos factores de expansión por distrito a partir de las estimaciones de población 2024 basadas en el censo 2017.

Trabajo preliminar

Paquetes R

library(sf)
library(sfnetworks)
library(dodgr)
library(tidygraph)
library(data.table)
library(dplyr)
library(purrr)
library(TSP)
library(spatstat)
library(openxlsx)
library(ggplot2)
library(ggspatial)
library(knitr)
library(tidyr)
library(stringr)
library(RColorBrewer)
library(foreach)
library(doFuture)
library(RColorBrewer)

Cargando los datos del cuestionario

Hogares <- read.xlsx("BD_Domiciliaria_Expandida_Procesada_HT.xlsx", sheet = "Hogares")
Personas <- read.xlsx("BD_Domiciliaria_Expandida_Procesada_HT.xlsx", sheet = "Personas")
viajes <- read.xlsx("BD_Domiciliaria_Expandida_Procesada_HT.xlsx", sheet = "Viajes")

EF <- read.xlsx("factores de emision lima 2024.xlsx")
Censo <- read.xlsx("Distritos_pop.xlsx")

Cargando los datos SHP

Nota: se sigue usando la red vial de 2010 que es la más desarrollada y extensa, sobre todo en zonas periféricas.

Zonas <- st_read(dsn = ".", layer = "zonas") %>% st_transform(4326)
## Reading layer `zonas' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima 2023' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1209 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.7979 ymin: -13.22341 xmax: -76.20938 ymax: -10.70688
## Geodetic CRS:  WGS 84
Zonas_urba <- st_read(dsn = ".", layer = "Zonas_urba") %>% st_transform(4326)
## Reading layer `Zonas_urba' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima 2023' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1185 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.19635 ymin: -12.4967 xmax: -76.66755 ymax: -11.72564
## Geodetic CRS:  WGS 84
setwd("C:/Users/hugot/Documents/Rennes 2/Calcul distance SIG/Reprise perso/RMD Lima")

manzanas <- st_read(dsn = "Data", layer = "LM_Manzanas_2020")
## Reading layer `LM_Manzanas_2020' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 102406 features and 27 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.19568 ymin: -12.49489 xmax: -76.66937 ymax: -11.72745
## Geodetic CRS:  WGS 84
limites_lima <- st_read(dsn = "Data", layer = "distritos") %>% st_transform(4326)
## Reading layer `distritos' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 49 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 254901.1 ymin: 8615341 xmax: 323611.5 ymax: 8719846
## Projected CRS: WGS 84 / UTM zone 18S
conos_lima <- st_read(dsn = "Data", layer = "Conos") %>% st_transform(4326)
## Reading layer `Conos' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.25154 ymin: -12.51954 xmax: -76.62155 ymax: -11.57301
## Geodetic CRS:  WGS 84
Routes <- st_read(dsn = "Data", layer = "Reseau_routier_2010")                            # Road network
## Reading layer `Reseau_routier_2010' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 252614 features and 20 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -77.19567 ymin: -12.49472 xmax: -76.50737 ymax: -11.72872
## Geodetic CRS:  WGS 84
Routes <- st_cast(Routes, "LINESTRING") %>% st_transform(4326)
Metro <- st_read(dsn = "Data", layer = "Metro1_connected")%>% st_transform(4326)          # Metro line 1
## Reading layer `Metro1_connected' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -77.01443 ymin: -12.20793 xmax: -76.93301 ymax: -11.95878
## Geodetic CRS:  WGS 84
Metro_stops <- st_read(dsn = "Data", layer = "Metro_estaciones")%>% st_transform(4326)    # Metro line 1 stops
## Reading layer `Metro_estaciones' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 26 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -77.01318 ymin: -12.20793 xmax: -76.93301 ymax: -11.95878
## Geodetic CRS:  WGS 84
BRT <- st_read(dsn = "Data", layer = "MetropolitanoAlimentador")%>% st_transform(4326)    # BRT network
## Reading layer `MetropolitanoAlimentador' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 132 features and 9 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -77.11262 ymin: -12.22851 xmax: -76.92755 ymax: -11.86455
## Geodetic CRS:  WGS 84
BRT_stops <- st_read(dsn = "Data", layer = "MetropolitanoAlimentador_estaciones")%>% st_transform(4326) #BRT stops 
## Reading layer `MetropolitanoAlimentador_estaciones' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 351 features and 16 fields
## Geometry type: MULTIPOINT
## Dimension:     XY
## Bounding box:  xmin: -77.11262 ymin: -12.22851 xmax: -76.92755 ymax: -11.86455
## Geodetic CRS:  WGS 84
Bus <- st_read(dsn = "Data", layer = "Bus")%>% st_transform(4326)                         # Bus lines
## Reading layer `Bus' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 142 features and 25 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -77.16914 ymin: -12.24639 xmax: -76.75636 ymax: -11.75389
## Geodetic CRS:  WGS 84
Custer <- st_read(dsn = "Data", layer = "Custer")%>% st_transform(4326)                   # Custer lines
## Reading layer `Custer' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 300 features and 25 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -77.17267 ymin: -12.49406 xmax: -76.66716 ymax: -11.75247
## Geodetic CRS:  WGS 84
Combi <- st_read(dsn = "Data", layer = "Combi") %>% st_transform(4326)                    # Combi lines
## Reading layer `Combi' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Lima\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 260 features and 25 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -77.14295 ymin: -12.38435 xmax: -76.77376 ymax: -11.72623
## Geodetic CRS:  WGS 84

Procesos preliminares

Se definen las categorías de modos que se usarán dentro del modelo de distancias. También se define una categoría para la comparación con los resultados para Bogotá.

modos <- data.frame(modo_principal_jer = unique(viajes$modo_principal_jer),
                    modo_principal_typo = c("Custer", "A_pie", "Taxi", "Corredor", "Mototaxi", "Moto", "Combi", "Patineta", "Auto", "Bus", "Alimentador", "Metropolitano", "Tren", "Otro", "Colectivo", "Bicicleta", "Taxi_por_aplicacion", "Trailer", "Camion pequeno", "Camion"),
                    modo_principal_comparado = c("Transporte público", "A pie", "Taxi / Carro por aplicación", "Transporte público", "Transporte informal", "Moto", "Transporte informal", "Otro", "Auto", "Transporte público", "Transporte público", "Transporte público", "Transporte público", "Otro", "Transporte informal", "Bicicleta", "Taxi / Carro por aplicación", "Otro", "Otro", "Otro"),
                    modo_principal_LATS = c("Regular bus", "Walking", "Taxi", "Regular bus", "Paratransit", "Moto", "Regular bus", "Other", "Private car", "Regular bus", "Regular bus", "BRT", "Rapid transit", "Other", "Paratransit", "Bike", "Taxi", "Other", "Other", "Other"))

viajes <- viajes %>% left_join(modos, by = "modo_principal_jer")

Excluimos los viajes entrando o saliendo del area conformado por la Provincia de Lima, para la cual tenemos un archivo SHP con la geometría de las zonas. A la fecha, se usa el archivo con las 1,209 Zonas de 2024, y se reduce el número de viajes en la base de 32,448 a 32,373 (-0.2%).

viajes <- viajes[viajes$DES_ZONA %in% Zonas$ID_ZONA & viajes$ORI_ZONA %in% Zonas$ID_ZONA,]

Necesitamos restringir la base de hogares a los 6 010 que efectivamente contestaron la encuesta, y luego crear un shapefile a partir de las coordenadas de su lugar de residencia.

viajes$duracion <- viajes$q24_1_hora_llegada+viajes$q24_2_minuto_llegada/60 - (viajes$q21_1_hora_salida+viajes$q21_2_minuto_salida/60)
# (quita 1 observacion)
viajes <- viajes[!is.na(viajes$DESLat) & !is.na(viajes$DESLon) & !is.na(viajes$ORILat) & !is.na(viajes$ORILon),]
Hogares <- Hogares[Hogares$filtro_ini == "Sí" & !(is.na(Hogares$filtro_ini)),]
Hogares$fexp_hogar <- as.numeric(Hogares$fexp_hogar)

Hogares <- Hogares %>% st_as_sf(coords = c("longitude", "latitude"), crs = 4326) 


#st_write(Hogares, "Hogares_loc.shp")

#Definición de la ventana
zoom_level <- 9.5
lon_span <- 360 / 2^zoom_level
lat_span <- 360 / 2^zoom_level

zoom_to_Lima <- c(-76.94, -12.05)  # Lima
lon_bounds_Lima <- c(zoom_to_Lima[1] - lon_span / 2, zoom_to_Lima[1] + lon_span / 2)
lat_bounds_Lima <- c(zoom_to_Lima[2] - lat_span / 2, zoom_to_Lima[2] + lat_span / 2)

ggplot()+
  theme_bw()+
  geom_sf(data = Zonas, fill = "white")+
  geom_sf(data = Hogares, col = "red", size = 0.5)+
  labs(title = "Ubicación de los hogares de la muestra", 
       subtitle = "Provincia de Lima", 
       caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, ATU, 2023", 
       col = "UPZ") +
  coord_sf(xlim = lon_bounds_Lima, ylim = lat_bounds_Lima, datum = NA)+
  labs(x = "", y = "") +
  theme(legend.position = "right",
        legend.text = element_text(size=10, hjust = 1),
        legend.title = element_text(size=15),
        plot.title = element_text(size=15),
        plot.caption = element_text(size = 7, face = "italic", hjust = 0, vjust = 12))+
  annotation_scale(location = "bl", height = unit(0.12, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(0.8, "cm"))

Finalmente necesitamos identificar para cada viaje si sus puntos de origen y destino son el hogar de la persona u otro. Para hacerlo, necesitamos buscar en la base de viajes las variables que contienen estas informaciones. Este procedimiento busca aumentar la precisión de los viajes que comienzan o terminan en el hogar.

#Lugar origen
colnames(viajes)[colnames(viajes) == "q19_lugar_origen"] <- "lugar_origen"
viajes$lugar_origen[viajes$lugar_origen != "Residencial (CASA / HOGAR)"] <- "Otro lugar"
viajes$lugar_origen[viajes$lugar_origen == "Residencial (CASA / HOGAR)"] <- "Hogar"

#Lugar destino
colnames(viajes)[colnames(viajes) == "q22_lugar_destino"] <- "lugar_destino"
viajes$lugar_destino[viajes$lugar_destino != "Residencial (CASA / HOGAR)"] <- "Otro lugar"
viajes$lugar_destino[viajes$lugar_destino == "Residencial (CASA / HOGAR)"] <- "Hogar"

#Control (no debe haber viajes Hogar --> Hogar)
viajes$test_origen_destino[viajes$lugar_origen == "Hogar" & viajes$lugar_destino == "Hogar"] <- "HOGAR_HOGAR"
viajes$test_origen_destino[viajes$lugar_origen == "Hogar" & viajes$lugar_destino == "Otro lugar"] <- "HOGAR_Otro"
viajes$test_origen_destino[viajes$lugar_origen == "Otro lugar" & viajes$lugar_destino == "Hogar"] <- "Otro_HOGAR"
viajes$test_origen_destino[viajes$lugar_origen == "Otro lugar" & viajes$lugar_destino == "Otro lugar"] <- "Otro_Otro"
print(table(viajes$test_origen_destino))
## 
## HOGAR_Otro Otro_HOGAR  Otro_Otro 
##      15666      15540       1166

Inicialización de las redes

Red vial

Routes$OBJECTID = NULL

# Ponderación de los segmentos y restricción al principal componente conexo
net_mot <- weight_streetnet(Routes, wt_profile = "motorcar", type_col = "type")
net_mot <- net_mot[net_mot$component == 1,]
net_cycl <- weight_streetnet(Routes, wt_profile = "bicycle", type_col = "type")
net_cycl <- net_cycl[net_cycl$component == 1,]
net_walk <- weight_streetnet(Routes, wt_profile = "foot", type_col = "type")
net_walk <- net_walk[net_walk$component == 1,]

Metro line 1

#Building the network
net_metro = as_sfnetwork(Metro, directed = FALSE) %>%
  st_transform(4326) %>%
  activate("edges") %>%
  mutate(weight = edge_length())

#Including (blend) the stations to the network
net_metro = st_network_blend(net_metro, Metro_stops)

#Computing the new length of each edge 
net_metro <- net_metro %>% activate("edges") %>% mutate(weight = edge_length())

BRT

BRT <- st_cast(BRT, "LINESTRING")
BRT_stops <- st_cast(BRT_stops, "POINT")

#Building the network
net_brt = as_sfnetwork(BRT, directed = FALSE) %>%
  st_transform(4326) %>%
  activate("edges") %>%
  mutate(weight = edge_length())

#Subdividing edges
net_brt = convert(net_brt, to_spatial_subdivision)

#Including (blend) the stations to the network
net_brt = st_network_blend(net_brt, BRT_stops)

net_brt <- net_brt %>%
  activate("nodes") %>%
  filter(group_components() == 1)

#Computing the new length of each edge 
net_brt <- net_brt %>% activate("edges") %>% mutate(weight = edge_length())

Bus

Bus <- st_cast(Bus, "LINESTRING")

#Building the network
net_bus = as_sfnetwork(Bus, directed = FALSE) %>%
  st_transform(4326) %>%
  activate("edges") %>%
  mutate(weight = edge_length())

#Subdividing edges
net_bus = convert(net_bus, to_spatial_subdivision)

#Computing the new length of each edge 
net_bus <- net_bus %>% activate("edges") %>% mutate(weight = edge_length())

Custer

Custer <- st_cast(Custer, "LINESTRING")

#Building the network
net_custer = as_sfnetwork(Custer, directed = FALSE) %>%
  st_transform(4326) %>%
  activate("edges") %>%
  mutate(weight = edge_length())

#Subdividing edges
net_custer = convert(net_custer, to_spatial_subdivision)

#Computing the new length of each edge
net_custer <- net_custer %>% activate("edges") %>% mutate(weight = edge_length())

Combi

Combi <- st_cast(Combi, "LINESTRING")

#Building the network
net_combi = as_sfnetwork(Combi, directed = FALSE) %>%
  st_transform(4326) %>%
  activate("edges") %>%
  mutate(weight = edge_length()) %>%
  activate("nodes") 

#Subdividing edges
net_combi = convert(net_combi, to_spatial_subdivision)

net_combi <- net_combi %>%
  activate("nodes") %>%
  filter(group_components() == 1)

#Computing the new length of each edge
net_combi <- net_combi %>% activate("edges") %>% mutate(weight = edge_length())

Zonificación

Creamos una matriz cuadrangular con malla fina (200 m) y agrupamos las coordenadas de los puntos de origen y destino a nivel de las zonas cuadradas. Este proceso busca acelerar el cálculo de las distancias respeto al metodo de cálculo de punto a punto.

#Proyección de la capa de zonas
Zonas_proj <- Zonas %>% st_transform(32718)

#creacion de capas geometricas con los puntos de partida y de llegada y conversion en coordenadas proyectadas
ori_viajes <- viajes[,c("id_viaje", "ORILon", "ORILat")] %>% st_as_sf(coords = c("ORILon", "ORILat"), crs = 4326) %>% st_transform(32718)
dest_viajes <- viajes[,c("id_viaje", "DESLon", "DESLat")] %>% st_as_sf(coords = c("DESLon", "DESLat"), crs = 4326) %>% st_transform(32718)
Malla (metros) Matriz cuadrangular (num puntos) Matriz reducida (num puntos)
5000 1960 92
500 194250 2647
200 1213625 6486 
#Creación de una matriz cuadrangular

cellsize <- 200

matriz <- st_make_grid(
  Zonas_proj,
  cellsize = cellsize,
  crs = 32718,
  what = "centers",
  square = TRUE)

matriz <- st_as_sf(matriz)
matriz$id_grid <- c(1:dim(matriz)[1])
id_ppv_ori <- data.frame(id_grid = st_nearest_feature(ori_viajes, matriz))
id_ppv_ori$id_viaje <- ori_viajes$id_viaje
id_ppv_dest <- data.frame(id_grid = st_nearest_feature(dest_viajes, matriz))
id_ppv_dest$id_viaje <- dest_viajes$id_viaje
matriz <- matriz[matriz$id_grid %in% id_ppv_ori$id_grid | matriz$id_grid %in% id_ppv_dest$id_grid,]

matriz <- matriz %>% st_transform(4326)
#viajes <- viajes[,1:110]
viajes <- viajes %>% left_join(id_ppv_ori, by = "id_viaje")
colnames(viajes)[colnames(viajes) == "id_grid"] <- "id_grid_ori"

viajes <- viajes %>% left_join(id_ppv_dest, by = "id_viaje")
colnames(viajes)[colnames(viajes) == "id_grid"] <- "id_grid_dest"

Cálculo de las distancias

Se calcula la longitud de cada viaje de acuerdo a los siguiente:

Modo principal Red empleada
Tren eléctrico Tren eléctrico
Metropolitano Metropolitano
Alimentador Metropolitano
Corredores Bus
Bus Bus
Custer Custer
Combi Combi
A pie Red vial (ponderación peatón)
Bicicleta Red vial (ponderación ciclista)
Cualquier otro Red vial (ponderación automovilista) 
print(Sys.time())

#5-6 min con malla 500m
#10 min con malla 200m

#LINEA RECTA (para sustituir NA)
mat_straight_22 <- as.data.frame(st_distance(x = st_geometry(matriz), y = st_geometry(matriz)))
rownames(mat_straight_22) <- matriz$id_grid
colnames(mat_straight_22) <- matriz$id_grid
mat_straight_22_table <- as.data.frame(as.table(as.matrix(mat_straight_22)))
colnames(mat_straight_22_table) <- c("id_grid_ori", "id_grid_dest", "dist_straight")

#TREN
mat_metro_22 <- st_network_cost(net_metro, from = st_geometry(matriz), to = st_geometry(matriz))
rownames(mat_metro_22) <- matriz$id_grid
colnames(mat_metro_22) <- matriz$id_grid
mat_metro_22_table <- as.data.frame(as.table(mat_metro_22))
colnames(mat_metro_22_table) <- c("id_grid_ori", "id_grid_dest", "dist_metro")

#BRT
mat_brt_22 <- st_network_cost(net_brt, from = st_geometry(matriz), to = st_geometry(matriz))
rownames(mat_brt_22) <- matriz$id_grid
colnames(mat_brt_22) <- matriz$id_grid
mat_brt_22_table <- as.data.frame(as.table(mat_brt_22))
colnames(mat_brt_22_table) <- c("id_grid_ori", "id_grid_dest", "dist_brt")

#BUS
mat_bus_22 <- st_network_cost(net_bus, from = st_geometry(matriz), to = st_geometry(matriz))
rownames(mat_bus_22) <- matriz$id_grid
colnames(mat_bus_22) <- matriz$id_grid
mat_bus_22_table <- as.data.frame(as.table(mat_bus_22))
colnames(mat_bus_22_table) <- c("id_grid_ori", "id_grid_dest", "dist_bus")

#CUSTER
mat_custer_22 <- st_network_cost(net_custer, from = st_geometry(matriz), to = st_geometry(matriz))
rownames(mat_custer_22) <- matriz$id_grid
colnames(mat_custer_22) <- matriz$id_grid
mat_custer_22_table <- as.data.frame(as.table(mat_custer_22))
colnames(mat_custer_22_table) <- c("id_grid_ori", "id_grid_dest", "dist_custer")

#COMBI
mat_combi_22 <- st_network_cost(net_combi, from = st_geometry(matriz), to = st_geometry(matriz))
rownames(mat_combi_22) <- matriz$id_grid
colnames(mat_combi_22) <- matriz$id_grid
mat_combi_22_table <- as.data.frame(as.table(mat_combi_22))
colnames(mat_combi_22_table) <- c("id_grid_ori", "id_grid_dest", "dist_combi")

#MOTORIZADO
mat_mot_22 <- dodgr_dists(graph = net_mot, from = as.data.frame(st_coordinates(matriz)), to = as.data.frame(st_coordinates(matriz)), shortest = TRUE)
rownames(mat_mot_22) <- matriz$id_grid
colnames(mat_mot_22) <- matriz$id_grid
mat_mot_22[is.na(mat_mot_22)] <- mat_straight_22[is.na(mat_mot_22)]
mat_mot_22_table <- as.data.frame(as.table(mat_mot_22))
colnames(mat_mot_22_table) <- c("id_grid_ori", "id_grid_dest", "dist_car")

#BICI
mat_cycl_22 <- dodgr_dists(graph = net_cycl, from = as.data.frame(st_coordinates(matriz)), to = as.data.frame(st_coordinates(matriz)), shortest = TRUE)
rownames(mat_cycl_22) <- matriz$id_grid
colnames(mat_cycl_22) <- matriz$id_grid
mat_cycl_22[is.na(mat_cycl_22)] <- mat_straight_22[is.na(mat_cycl_22)]
mat_cycl_22_table <- as.data.frame(as.table(mat_cycl_22))
colnames(mat_cycl_22_table) <- c("id_grid_ori", "id_grid_dest", "dist_bike")

#WALK
mat_walk_22 <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(matriz)), to = as.data.frame(st_coordinates(matriz)), shortest = TRUE)
rownames(mat_walk_22) <- matriz$id_grid
colnames(mat_walk_22) <- matriz$id_grid
mat_walk_22[is.na(mat_walk_22)] <- mat_straight_22[is.na(mat_walk_22)]
mat_walk_22_table <- as.data.frame(as.table(mat_walk_22))
colnames(mat_walk_22_table) <- c("id_grid_ori", "id_grid_dest", "dist_walk")

print(Sys.time())
viajes_dist <- viajes[,c("id_viaje", "id_per", "id_hogar", "lugar_origen", "lugar_destino", "id_grid_ori", "id_grid_dest", "fexp_viaje", "duracion", "modo_principal_typo", "modo_principal_comparado", "test_origen_destino")]

#Unión por columnas múltiples 
#1 min con malla 500m
#6 min con malla 200m

print(Sys.time())

viajes_dist <- merge(viajes_dist, mat_mot_22_table, by = c("id_grid_ori", "id_grid_dest")) 

viajes_dist <- merge(viajes_dist, mat_cycl_22_table, by = c("id_grid_ori", "id_grid_dest")) 

viajes_dist <- merge(viajes_dist, mat_walk_22_table, by = c("id_grid_ori", "id_grid_dest")) 

viajes_dist <- merge(viajes_dist, mat_metro_22_table, by = c("id_grid_ori", "id_grid_dest"))

viajes_dist <- merge(viajes_dist, mat_brt_22_table, by = c("id_grid_ori", "id_grid_dest")) 

viajes_dist <- merge(viajes_dist, mat_bus_22_table, by = c("id_grid_ori", "id_grid_dest")) 

viajes_dist <- merge(viajes_dist, mat_custer_22_table, by = c("id_grid_ori", "id_grid_dest")) 

viajes_dist <- merge(viajes_dist, mat_combi_22_table, by = c("id_grid_ori", "id_grid_dest")) 

viajes_dist <- merge(viajes_dist, mat_straight_22_table, by = c("id_grid_ori", "id_grid_dest"))

print(Sys.time())

write.xlsx(viajes_dist, "viajes_dist_no_corr_matriz_200.xlsx")
viajes_dist <- read.xlsx("Tablas producidas/viajes_dist_no_corr_matriz_200.xlsx")

viajes_dist$Distancia <- viajes_dist$dist_car

viajes_dist$Distancia[viajes_dist$modo_principal_typo == "A_pie"] <- viajes_dist$dist_walk[viajes_dist$modo_principal_typo == "A_pie"]

viajes_dist$Distancia[viajes_dist$modo_principal_typo %in% c("Bicicleta")] <- viajes_dist$dist_bike[viajes_dist$modo_principal_typo %in% c("Bicicleta")]

viajes_dist$Distancia[viajes_dist$modo_principal_typo %in% c("Alimentador", "Metropolitano")] <- viajes_dist$dist_brt[viajes_dist$modo_principal_typo %in% c("Alimentador", "Metropolitano")]

viajes_dist$Distancia[viajes_dist$modo_principal_typo == "Tren"] <- viajes_dist$dist_metro[viajes_dist$modo_principal_typo == "Tren"]

viajes_dist$Distancia[viajes_dist$modo_principal_typo %in% c("Bus", "Corredor")] <- viajes_dist$dist_bus[viajes_dist$modo_principal_typo %in% c("Bus", "Corredor")]

viajes_dist$Distancia[viajes_dist$modo_principal_typo == "Custer"] <- viajes_dist$dist_custer[viajes_dist$modo_principal_typo == "Custer"]

viajes_dist$Distancia[viajes_dist$modo_principal_typo == "Combi"] <- viajes_dist$dist_combi[viajes_dist$modo_principal_typo == "Combi"]

viajes_dist$Distancia[is.na(viajes_dist$Distancia)] <- viajes_dist$dist_straight[is.na(viajes_dist$Distancia)]

viajes_dist$Distancia <- as.numeric(viajes_dist$Distancia)

Para los viajes cuyo modo principal es el transporte colectivo (Tren, BRT, bus, combi, cúster), añadimos las distancias hasta y desde los paraderos.

plot(net_metro)

#Conversion de nudos de la red en formato sf. Incluye las 26 estaciones
acceso_metro <- net_metro %>%
  activate("nodes") %>%
  st_geometry() %>%
  st_as_sf()

#Matriz de distancias a pie entre puntos de la matriz y estaciones
mat_walk_matriz_metro <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(matriz)), to = as.data.frame(st_coordinates(acceso_metro)), shortest = TRUE)
rownames(mat_walk_matriz_metro) <- matriz$id_grid
colnames(mat_walk_matriz_metro) <- rownames(acceso_metro)

#Búsqueda de la estación de metro más cercana a cada punto de la matriz
distancia <- apply(mat_walk_matriz_metro[,1:dim(mat_walk_matriz_metro)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_matriz_metro)[apply(mat_walk_matriz_metro, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_matriz <- cbind(name, distancia)
estacion_mas_cerca_matriz$id_grid <- as.numeric(rownames(estacion_mas_cerca_matriz))

#Adición de la distancia hasta y desde el metro en los viajes con este modo
matriz_metro <- viajes_dist %>% left_join(estacion_mas_cerca_matriz[,c("id_grid", "distancia")], by = c("id_grid_ori" = "id_grid")) %>%
  group_by(id_grid_ori) %>%
  summarize(dist_matriz_metro = unique(distancia))

metro_matriz <- viajes_dist %>% left_join(estacion_mas_cerca_matriz[,c("id_grid", "distancia")], by = c("id_grid_dest" = "id_grid")) %>%
  group_by(id_grid_dest) %>%
  summarize(dist_metro_matriz = unique(distancia))

viajes_dist <- viajes_dist %>% 
  left_join(matriz_metro, by = "id_grid_ori") %>%
  left_join(metro_matriz, by = "id_grid_dest")

viajes_dist$dist_metro_incl_caminata <- as.numeric(viajes_dist$dist_metro) + viajes_dist$dist_matriz_metro + viajes_dist$dist_metro_matriz
plot(net_brt)

#Conversion de nudos de la red en formato sf. Incluye las 351 estaciones y paraderos del Metropolitano+Alimentador más 1 070 nudos de las rutas de Alimentación que no pudieron ser evitados en la construcción de la red.
acceso_brt <- net_brt %>%
  activate("nodes") %>%
  st_geometry() %>%
  st_as_sf()

#Matriz de distancias a pie entre puntos de la matriz y estaciones
mat_walk_matriz_brt <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(matriz)), to = as.data.frame(st_coordinates(acceso_brt)), shortest = TRUE)
rownames(mat_walk_matriz_brt) <- matriz$id_grid
colnames(mat_walk_matriz_brt) <- rownames(acceso_brt)

#Búsqueda de la estación de brt más cercana a cada punto de la matriz
distancia <- apply(mat_walk_matriz_brt[,1:dim(mat_walk_matriz_brt)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_matriz_brt)[apply(mat_walk_matriz_brt, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_matriz <- cbind(name, distancia)
estacion_mas_cerca_matriz$id_grid <- as.numeric(rownames(estacion_mas_cerca_matriz))

#Adición de la distancia hasta y desde el brt en los viajes con este modo
matriz_brt <- viajes_dist %>% left_join(estacion_mas_cerca_matriz[,c("id_grid", "distancia")], by = c("id_grid_ori" = "id_grid")) %>%
  group_by(id_grid_ori) %>%
  summarize(dist_matriz_brt = unique(distancia))

brt_matriz <- viajes_dist %>% left_join(estacion_mas_cerca_matriz[,c("id_grid", "distancia")], by = c("id_grid_dest" = "id_grid")) %>%
  group_by(id_grid_dest) %>%
  summarize(dist_brt_matriz = unique(distancia))

viajes_dist <- viajes_dist %>% 
  left_join(matriz_brt, by = "id_grid_ori") %>%
  left_join(brt_matriz, by = "id_grid_dest")

viajes_dist$dist_brt_incl_caminata <- as.numeric(viajes_dist$dist_brt) + viajes_dist$dist_matriz_brt + viajes_dist$dist_brt_matriz
#Conversion de nudos de la red en formato sf. Incluye los 12 192 nudos de la red de bus.
acceso_bus <- net_bus %>%
  activate("nodes") %>%
  st_geometry() %>%
  st_as_sf()

plot(acceso_bus)

#Matriz de distancias a pie entre puntos de la matriz y estaciones
mat_walk_matriz_bus <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(matriz)), to = as.data.frame(st_coordinates(acceso_bus)), shortest = TRUE)
rownames(mat_walk_matriz_bus) <- matriz$id_grid
colnames(mat_walk_matriz_bus) <- rownames(acceso_bus)

#Búsqueda de la estación de bus más cercana a cada punto de la matriz
distancia <- apply(mat_walk_matriz_bus[,1:dim(mat_walk_matriz_bus)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_matriz_bus)[apply(mat_walk_matriz_bus, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_matriz <- cbind(name, distancia)
estacion_mas_cerca_matriz$id_grid <- as.numeric(rownames(estacion_mas_cerca_matriz))

#Adición de la distancia hasta y desde el bus en los viajes con este modo
matriz_bus <- viajes_dist %>% left_join(estacion_mas_cerca_matriz[,c("id_grid", "distancia")], by = c("id_grid_ori" = "id_grid")) %>%
  group_by(id_grid_ori) %>%
  summarize(dist_matriz_bus = unique(distancia))

bus_matriz <- viajes_dist %>% left_join(estacion_mas_cerca_matriz[,c("id_grid", "distancia")], by = c("id_grid_dest" = "id_grid")) %>%
  group_by(id_grid_dest) %>%
  summarize(dist_bus_matriz = unique(distancia))

viajes_dist <- viajes_dist %>% 
  left_join(matriz_bus, by = "id_grid_ori") %>%
  left_join(bus_matriz, by = "id_grid_dest")

viajes_dist$dist_bus_incl_caminata <- as.numeric(viajes_dist$dist_bus) + viajes_dist$dist_matriz_bus + viajes_dist$dist_bus_matriz
#Conversion de nudos de la red en formato sf. Incluye los 15 989 nudos de la red de custer.
acceso_custer <- net_custer %>%
  activate("nodes") %>%
  st_geometry() %>%
  st_as_sf()

plot(acceso_custer)

#Matriz de distancias a pie entre puntos de la matriz y estaciones
mat_walk_matriz_custer <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(matriz)), to = as.data.frame(st_coordinates(acceso_custer)), shortest = TRUE)
rownames(mat_walk_matriz_custer) <- matriz$id_grid
colnames(mat_walk_matriz_custer) <- rownames(acceso_custer)

#Búsqueda de la estación de custer más cercana a cada punto de la matriz
distancia <- apply(mat_walk_matriz_custer[,1:dim(mat_walk_matriz_custer)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_matriz_custer)[apply(mat_walk_matriz_custer, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_matriz <- cbind(name, distancia)
estacion_mas_cerca_matriz$id_grid <- as.numeric(rownames(estacion_mas_cerca_matriz))

#Adición de la distancia hasta y desde el custer en los viajes con este modo
matriz_custer <- viajes_dist %>% left_join(estacion_mas_cerca_matriz[,c("id_grid", "distancia")], by = c("id_grid_ori" = "id_grid")) %>%
  group_by(id_grid_ori) %>%
  summarize(dist_matriz_custer = unique(distancia))

custer_matriz <- viajes_dist %>% left_join(estacion_mas_cerca_matriz[,c("id_grid", "distancia")], by = c("id_grid_dest" = "id_grid")) %>%
  group_by(id_grid_dest) %>%
  summarize(dist_custer_matriz = unique(distancia))

viajes_dist <- viajes_dist %>% 
  left_join(matriz_custer, by = "id_grid_ori") %>%
  left_join(custer_matriz, by = "id_grid_dest")

viajes_dist$dist_custer_incl_caminata <- as.numeric(viajes_dist$dist_custer) + viajes_dist$dist_matriz_custer + viajes_dist$dist_custer_matriz
#Conversion de nudos de la red en formato sf. Incluye los 11 797 nudos de la red de combi.
acceso_combi <- net_combi %>%
  activate("nodes") %>%
  st_geometry() %>%
  st_as_sf()

plot(acceso_combi)

#Matriz de distancias a pie entre puntos de la matriz y estaciones
mat_walk_matriz_combi <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(matriz)), to = as.data.frame(st_coordinates(acceso_combi)), shortest = TRUE)
rownames(mat_walk_matriz_combi) <- matriz$id_grid
colnames(mat_walk_matriz_combi) <- rownames(acceso_combi)

#Búsqueda de la estación de combi más cercana a cada punto de la matriz
distancia <- apply(mat_walk_matriz_combi[,1:dim(mat_walk_matriz_combi)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_matriz_combi)[apply(mat_walk_matriz_combi, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_matriz <- cbind(name, distancia)
estacion_mas_cerca_matriz$id_grid <- as.numeric(rownames(estacion_mas_cerca_matriz))

#Adición de la distancia hasta y desde el combi en los viajes con este modo
matriz_combi <- viajes_dist %>% left_join(estacion_mas_cerca_matriz[,c("id_grid", "distancia")], by = c("id_grid_ori" = "id_grid")) %>%
  group_by(id_grid_ori) %>%
  summarize(dist_matriz_combi = unique(distancia))

combi_matriz <- viajes_dist %>% left_join(estacion_mas_cerca_matriz[,c("id_grid", "distancia")], by = c("id_grid_dest" = "id_grid")) %>%
  group_by(id_grid_dest) %>%
  summarize(dist_combi_matriz = unique(distancia))

viajes_dist <- viajes_dist %>% 
  left_join(matriz_combi, by = "id_grid_ori") %>%
  left_join(combi_matriz, by = "id_grid_dest")

viajes_dist$dist_combi_incl_caminata <- as.numeric(viajes_dist$dist_combi) + viajes_dist$dist_matriz_combi + viajes_dist$dist_combi_matriz
viajes_dist$Distancia_incl_caminata <- viajes_dist$Distancia

viajes_dist$Distancia_incl_caminata[viajes_dist$modo_principal_typo == "Tren"] <- viajes_dist$dist_metro_incl_caminata[viajes_dist$modo_principal_typo == "Tren"]

viajes_dist$Distancia_incl_caminata[viajes_dist$modo_principal_typo %in% c("Alimentador", "Metropolitano")] <- viajes_dist$dist_brt_incl_caminata[viajes_dist$modo_principal_typo %in% c("Alimentador", "Metropolitano")]

viajes_dist$Distancia_incl_caminata[viajes_dist$modo_principal_typo %in% c("Bus", "Corredor")] <- viajes_dist$dist_bus_incl_caminata[viajes_dist$modo_principal_typo %in% c("Bus", "Corredor")]

viajes_dist$Distancia_incl_caminata[viajes_dist$modo_principal_typo == "Custer"] <- viajes_dist$dist_custer_incl_caminata[viajes_dist$modo_principal_typo == "Custer"]

viajes_dist$Distancia_incl_caminata[viajes_dist$modo_principal_typo == "Combi"] <- viajes_dist$dist_combi_incl_caminata[viajes_dist$modo_principal_typo == "Combi"]

Hallamos la distancia promedio desde cada red hasta los puntos de la matriz (con resolución 200m).

comparacion <- data.frame(red = c("Tren", "BRT", "Bus", "Custer", "Combi"),
                          dist_promedio_desde_matriz = c(mean(metro_matriz$dist_metro_matriz),
                                                         mean(brt_matriz$dist_brt_matriz),
                                                         mean(bus_matriz$dist_bus_matriz),
                                                         mean(custer_matriz$dist_custer_matriz),
                                                         mean(combi_matriz$dist_combi_matriz)),
                          dist_mediana_desde_matriz = c(median(metro_matriz$dist_metro_matriz),
                                                        median(brt_matriz$dist_brt_matriz),
                                                        median(bus_matriz$dist_bus_matriz),
                                                        median(custer_matriz$dist_custer_matriz),
                                                        median(combi_matriz$dist_combi_matriz)))

colnames(comparacion) <- c("Red", "Promedio", "Mediana")

kable(comparacion, caption = "Comparación de las distancias entre los puntos de la matriz (200m) y las redes de transporte colectivo", format.args = list(big.mark = ","))
Comparación de las distancias entre los puntos de la matriz (200m) y las redes de transporte colectivo
Red Promedio Mediana
Tren 12,829.877 9,080.7645
BRT 7,484.083 3,834.8841
Bus 1,905.864 301.2881
Custer 500.474 224.4569
Combi 2,463.278 336.1882

Para los viajes intrazona, se decide no abordarlos, dado que se requeriría tener un shape con las geometrías de las zonas con clase POLYGON y no POINT. Consideramos que la malla que escogimos es lo suficiente fina para que no haya un número importante de viajes intrazona. De hecho, representan 893 viajes (el 2.7% total).

Identificamos los viajes de ida y vuelta. Un viaje de ida y vuelta se identifica porque cumple las siguientes características:

Nota: dado que no pudimos identificar todos los viajes que inician en el hogar, sino solo el primero del día, este algoritmo solo puede detectar máximo una ida y vuelta por persona.

#Esta etapa es necesaria en el marco del control de velocidad de la siguiente etapa. Se define como ida y vuelta un trayecto casa -> cierto lugar -> casa realizado con el mismo modo de transporte. El lugar de destino del primer viaje y el lugar de origen del segundo deben coincidir. Si el modo de transporte cambia, no lo definimos como ida y vuelta.

ida <- viajes_dist[viajes_dist$lugar_origen == "Hogar", c("id_per", "id_viaje", "lugar_destino", "modo_principal_typo")]
vuelta <- viajes_dist[viajes_dist$lugar_destino == "Hogar", c("id_per", "id_viaje", "lugar_origen", "modo_principal_typo")]
ida_vuelta <- na.omit(ida %>% left_join(vuelta, by = "id_per"))
ida_vuelta$is_ida_vuelta <- (ida_vuelta$lugar_destino == ida_vuelta$lugar_origen) & (ida_vuelta$modo_principal_typo.x == ida_vuelta$modo_principal_typo.y)
ida <- ida_vuelta[,c("id_viaje.x", "is_ida_vuelta")]
colnames(ida)[colnames(ida) == "id_viaje.x"] <- "id_viaje"
ida <- ida %>% group_by(id_viaje) %>% summarize(is_ida_vuelta = as.logical(max(is_ida_vuelta)))
vuelta <- ida_vuelta[,c("id_viaje.y", "is_ida_vuelta")]
colnames(vuelta)[colnames(vuelta) == "id_viaje.y"] <- "id_viaje"
vuelta <- vuelta %>% group_by(id_viaje) %>% summarize(is_ida_vuelta = as.logical(max(is_ida_vuelta)))
viajes_dist <- viajes_dist %>% left_join(ida, by = c("id_viaje"))
viajes_dist <- viajes_dist %>% left_join(vuelta, by = c("id_viaje"))
viajes_dist$is_ida_vuelta.x[is.na(viajes_dist$is_ida_vuelta.x)] <- FALSE
viajes_dist$is_ida_vuelta.y[is.na(viajes_dist$is_ida_vuelta.y)] <- FALSE
viajes_dist$is_ida_vuelta <- viajes_dist$is_ida_vuelta.x | viajes_dist$is_ida_vuelta.y
viajes_dist$is_ida_vuelta.x <- NULL
viajes_dist$is_ida_vuelta.y <- NULL
viajes_dist$duracion_min <- 60*viajes_dist$duracion

#velocidad en km/h
viajes_dist$velocidad <- 60*viajes_dist$Distancia/(1000*viajes_dist$duracion_min)

#825 viajes a pie son demasiado rápidos (9%)
viajes_dist$Distancia[viajes_dist$modo_principal_typo == "A_pie" & viajes_dist$velocidad > 10] <- 1000*viajes_dist$duracion_min[viajes_dist$modo_principal_typo == "A_pie" & viajes_dist$velocidad > 10]*10/60

#5 viajes en bici son demasiado rápidos (2.7%)
viajes_dist$Distancia[viajes_dist$modo_principal_typo == "Bicicleta" & viajes_dist$velocidad > 30] <- 1000*viajes_dist$duracion_min[viajes_dist$modo_principal_typo == "Bicicleta" & viajes_dist$velocidad > 30]*30/60

#501 viajes en modos motorizados son demasiado rápidos (2.2%)
viajes_dist$Distancia[!(viajes_dist$modo_principal_typo %in% c("A_pie","Bicicleta")) & viajes_dist$velocidad > 50] <- 1000*viajes_dist$duracion_min[!(viajes_dist$modo_principal_typo %in% c("A_pie","Bicicleta")) & viajes_dist$velocidad > 50]*50/60

#se aplica una ultima correccion: los viajes de ida y vuelta deben tener la misma distancia. Se corrige la distancia de los viajes demasiado rápidos, pero en el caso de los viajes de ida y vuelta, dado que esperamos que la distancia coincida en la ida y el regreso, le vamos a asignar el maximo entre los dos valores.
ida_vuelta <- viajes_dist[viajes_dist$is_ida_vuelta == TRUE,] %>% 
  group_by(id_per, modo_principal_typo) %>%
  summarise(Distancia_AR = max(Distancia))
viajes_dist <- viajes_dist %>% left_join(ida_vuelta, by = c("id_per", "modo_principal_typo"))
viajes_dist$Distancia[viajes_dist$is_ida_vuelta == TRUE] <- viajes_dist$Distancia_AR[viajes_dist$is_ida_vuelta == TRUE]
viajes_dist$Distancia_AR <- NULL
viajes_dist$velocidad <- NULL

Por fin podemos crear la base final con las variables correspondientes a distancia, emisiones de gases de efecto invernadero y contaminantes aéreos.

viajes_dist <- viajes_dist %>% left_join(EF, by = "modo_principal_typo")
viajes_dist$v_moy = NULL

viajes_dist$`CO2-eq` <- viajes_dist$`CO2-eq`*viajes_dist$Distancia/1000
viajes_dist$CO <- viajes_dist$CO*viajes_dist$Distancia/1000
viajes_dist$NOx <- viajes_dist$NOx*viajes_dist$Distancia/1000
viajes_dist$BC <- viajes_dist$BC*viajes_dist$Distancia/1000
viajes_dist$COV <- viajes_dist$COV*viajes_dist$Distancia/1000
viajes_dist$PM <- viajes_dist$PM*viajes_dist$Distancia/1000

Actualización población

Evidenciamos problemas con los factores de expansión de la población de la encuesta 2023. Algunos distritos tienen sobrestimada su población total mientras que otros la tienen subestimada, respeto a la proyección del INEI para 2024.

También hay un problema con el muestreo, la mitad de las zonas no habiendo sido muestreadas. Como consecuencia, la encuesta no es respresentativa a nivel de zonas y solo puede alcanzar a serlo a nivel de distritos. Por lo tanto, adoptamos dos soluciones:

Hogares$distrito[Hogares$distrito == "BREÑA"] <- "BRENA"
Hogares$distrito[Hogares$distrito == "ATE"] <- "ATE VITARTE"
Hogares$distrito[Hogares$distrito == "MI PERÚ"] <- "MI PERU"

Hogares <- Hogares %>% st_join(Zonas[,c("ID_ZONA")])

Personas_indic <- Personas %>% 
  left_join(Hogares[,c("id_hog", "ID_ZONA")], by = "id_hog")%>%
  group_by(ID_ZONA) %>%
  summarise(pop = sum(fexp_per))

viajes_dist <- viajes_dist %>% left_join(Hogares[,c("id_hog", "distrito")], by = c("id_hogar" = "id_hog")) 

Zonas_urba <- Zonas_urba %>% left_join(Personas_indic, by = "ID_ZONA")
#calculamos la población por distrito segun el censo del 2017. Exportamos el archivo en formato Excel. Fuera de R, creamos una nueva columna con la población de los distritos proyectada al 2024.
Censo <- manzanas %>% st_drop_geometry() %>% group_by(DISTRITO) %>% summarize(pop_censo_2017 = sum(POBLACION))

#write.xlsx(Censo, "Distritos_pop.xlsx")
Censo <- read.xlsx("Distritos_pop.xlsx")

Comparación de la población por distrito segun segun encuesta 2023 y proyecciones INEI 2024

Distritos <- Zonas_urba %>% 
  group_by(distrito) %>% 
  summarise(pop_enc_2023 = sum(pop, na.rm = TRUE))

Distritos$distrito[Distritos$distrito == "BREŃA"] <- "BRENA"
Distritos$distrito[Distritos$distrito == "CARMEN DE LA LEGUA"] <- "CARMEN DE LA LEGUA REYNOSO"
Distritos <- Distritos %>% left_join(Censo[,c("DISTRITO", "pop_est_2024")], by = c("distrito" = "DISTRITO"))
Distritos$variacion <- 100*(Distritos$pop_enc_2023-Distritos$pop_est_2024)/Distritos$pop_est_2024

Distritos <- Distritos[order(Distritos$variacion, decreasing = TRUE),]

total <- data.frame("TOTAL", 
                   sum(Distritos$pop_enc_2023), 
                   sum(Distritos$pop_est_2024),
                   100*(sum(Distritos$pop_enc_2023)-sum(Distritos$pop_est_2024))/sum(Distritos$pop_est_2024))

Distritos_table <- st_drop_geometry(Distritos)

colnames(total) <- colnames(Distritos_table)

Distritos_table <- rbind(Distritos_table, total)

Distritos_table$pop_enc_2023 <- round(Distritos_table$pop_enc_2023, 0)

colnames(Distritos_table) <- c("Distrito", "Población encuesta ATU 2023", "Población proyectada INEI 2024", "% Variación")

kable(Distritos_table, caption = "Población por distrito", format.args = list(big.mark = ","))
Población por distrito
Distrito Población encuesta ATU 2023 Población proyectada INEI 2024 % Variación
SAN BARTOLO 147,378 8,600 1,613.6995035
SANTA MARIA DEL MAR 13,797 1,200 1,049.7459083
PUNTA NEGRA 67,910 8,300 718.1961590
PUCUSANA 124,918 17,300 622.0705694
LA PUNTA 32,201 4,600 600.0200717
SANTA ROSA 211,560 32,400 552.9637315
CIENEGUILLA 261,525 40,100 552.1812579
PUNTA HERMOSA 114,327 18,600 514.6619048
BARRANCO 159,587 41,200 287.3480177
CHACLACAYO 170,405 50,600 236.7687449
ANCON 224,168 73,600 204.5754319
SAN ISIDRO 220,218 73,100 201.2551484
SURQUILLO 270,781 108,400 149.7983755
CARMEN DE LA LEGUA REYNOSO 119,901 49,600 141.7362006
LURIN 245,137 102,200 139.8596937
PACHACAMAC 289,003 126,100 129.1855336
RIMAC 452,892 206,300 119.5310179
LA PERLA 156,370 73,000 114.2057518
LA MOLINA 355,705 167,400 112.4882064
MI PERU 98,706 52,400 88.3693376
BELLAVISTA 140,026 88,800 57.6872252
SAN BORJA 209,404 135,800 54.2005065
PUEBLO LIBRE 147,128 99,500 47.8674661
LURIGANCHO 372,739 280,800 32.7419289
EL AGUSTINO 276,789 233,000 18.7937128
LINCE 77,667 65,400 18.7575621
MIRAFLORES 134,325 116,500 15.3004580
CHORRILLOS 413,850 369,500 12.0027924
INDEPENDENCIA 256,145 248,300 3.1595559
SAN LUIS 62,321 61,500 1.3354899
CARABAYLLO 380,245 383,100 -0.7451670
VENTANILLA 355,133 365,600 -2.8630104
PUENTE PIEDRA 368,375 383,100 -3.8436644
LA VICTORIA 195,266 205,000 -4.7483249
SANTA ANITA 218,317 229,700 -4.9555242
SAN JUAN DE MIRAFLORES 389,919 418,000 -6.7180486
SAN MIGUEL 163,162 184,800 -11.7088453
VILLA MARIA DEL TRIUNFO 375,220 466,600 -19.5841869
MAGDALENA DEL MAR 56,416 72,000 -21.6444114
JESUS MARIA 67,666 90,200 -24.9824963
SANTIAGO DE SURCO 289,871 391,900 -26.0345509
LIMA 233,250 318,500 -26.7659861
BRENA 72,341 101,300 -28.5873096
VILLA EL SALVADOR 306,342 460,400 -33.4617291
CALLAO 336,568 535,300 -37.1253791
LOS OLIVOS 227,745 384,000 -40.6913311
ATE VITARTE 409,863 699,500 -41.4063186
SAN MARTIN DE PORRES 359,639 770,100 -53.2996481
COMAS 283,912 612,400 -53.6394429
SAN JUAN DE LURIGANCHO 337,287 1,216,300 -72.2693873
TOTAL 11,253,424 11,241,900 0.1025095 
#st_write(Distritos, "Distritos.shp")

Nuevos factores de expansión

Para actualizar el factor de expansión de las personas:

  • Le atribuimos un distrito de residencia a cada persona por unión con la base de hogares.
  • Hacemos la suma del número de personas entrevistadas dentro de cada distrito.
  • Finalmente, calculamos el factor de expansión de cada persona como el cociente entre el número de personas estimado dentro del distrito con el número de personas entrevistadas en el distrito.

Finalmente, todas las personas encuestadas dentro de un mismo distrito tendrán el mismo factor de expansión.

Para actualizar el factor de expansión de los hogares (no lo usaremos, solo será para fines informativos), le atribuimos a cada hogar el nuevo factor de expansión único de sus miembros. La nueva suma de los factores de expansión de los hogares ahora asciende a 3,570,012, mientras que la proyección del INEI para el 2024 es inferior con un valor de 2,813,100. Esto se debe a la calidad del muestreo inicial, que no pudo controlar a la vez la distribución de los hogares y la distribución de las personas dentro de los mismos.

Para actualizar el factor de expansión de los viajes, sustituimos el factor bruto por el nuevo factor de expansión de la persona y lo multiplicamos por el factor de ajuste modal.

#Distrito de referencia y cálculo del número de personas entrevistadas dentro de cada distrito
Personas <- Personas %>% left_join(Hogares[,c("id_hog", "distrito")], by = "id_hog")
Personas_por_distrito <- Personas %>% group_by(distrito) %>% summarise(num_pers_distrito = n()) %>%
  left_join(st_drop_geometry(Distritos[,c("distrito", "pop_est_2024")]), by = "distrito")

#Cálculo de los nuevos factores de expansión
Personas_por_distrito$fexp_per_nuevo <- Personas_por_distrito$pop_est_2024/Personas_por_distrito$num_pers_distrito
Personas <- Personas %>% left_join(Personas_por_distrito[,c("distrito", "fexp_per_nuevo")], by = "distrito")
#Cálculo del número de personas por hogar
Personas_por_hog <- Personas %>% group_by(id_hog) %>% summarise(num_pers_hog = n(), fexp_hogar_nuevo = unique(fexp_per_nuevo))
Hogares <- Hogares %>% left_join(Personas_por_hog, by = "id_hog")
#Unión con la base de personas para tener el factor de viaje bruto
viajes <- viajes %>% left_join(Personas[,c("id_pers", "fexp_per_nuevo")], by = c("id_per" = "id_pers"))

#Ajuste segun el modo
viajes$fexp_viaje_nuevo <- viajes$fexp_per_nuevo*viajes$f_ajust

viajes_dist <- viajes_dist[,1:47] %>% left_join(viajes[,c("id_viaje", "fexp_viaje_nuevo")], by = "id_viaje")

#write.xlsx(viajes_dist, paste0("viajes_dist_final_matriz_con_caminata_previa.xlsx"), overwrite = TRUE)

Estadísticas

Las estadísticas quedan así, con modos desagrupados y agrupados para comparación con Bogotá.

NOTA IMPORTANTE PARA LOS VIAJES CUYO MODO PRINCIPAL ES UN MODO COLECTIVO:

viajes_dist <- read.xlsx("Tablas producidas/viajes_dist_final_matriz_con_caminata_previa.xlsx")


estadisticas <- viajes_dist %>% 
  group_by(modo_principal_typo) %>% 
  summarize(numero_viajes = round(sum(fexp_viaje_nuevo),0), 
            mpkt = round(sum(fexp_viaje_nuevo*Distancia/1000000000),1),
            mpkt_incl_caminata = round(sum(fexp_viaje_nuevo*Distancia_incl_caminata/1000000000),1),
            reparto_modal = round(100*sum(fexp_viaje_nuevo)/sum(viajes_dist$fexp_viaje_nuevo),1), 
            reparto_modal_pkt = round(100*sum(fexp_viaje_nuevo*Distancia_incl_caminata)/sum(viajes_dist$fexp_viaje_nuevo*viajes_dist$Distancia_incl_caminata),1), 
            `CO2-eq`= round(sum(fexp_viaje_nuevo*`CO2-eq`/1000000),1),
            dist_por_viaj = round(sum(fexp_viaje_nuevo*Distancia_incl_caminata/1000)/sum(fexp_viaje_nuevo),1))
  
estadisticas <- estadisticas[order(estadisticas$numero_viajes, decreasing = TRUE),]

total <- data.frame("TOTAL", 
                   sum(estadisticas$numero_viajes),
                   sum(estadisticas$mpkt),
                   sum(estadisticas$mpkt_incl_caminata),
                   sum(estadisticas$reparto_modal),
                   sum(estadisticas$reparto_modal_pkt),
                   sum(estadisticas$`CO2-eq`),
                   round(sum(estadisticas$numero_viajes*estadisticas$dist_por_viaj)/sum(estadisticas$numero_viajes),1))

colnames(total) <- colnames(estadisticas)

estadisticas <- rbind(estadisticas, total)

colnames(estadisticas) <- c("Modo principal","Número de viajes","Millón PKT (solo tramo principal)", "Millón PKT (total incluyendo caminata inicial y final)","% Número de viajes","% PKT","Emisiones GEI (tCO2-eq)", "Distancia media por viaje (km)")

kable(estadisticas, caption = "Estadísticas básicas", format.args = list(big.mark = ","))
Estadísticas básicas
Modo principal Número de viajes Millón PKT (solo tramo principal) Millón PKT (total incluyendo caminata inicial y final) % Número de viajes % PKT Emisiones GEI (tCO2-eq) Distancia media por viaje (km)
A_pie 5,500,873 4.8 4.8 22.7 2.4 0.0 0.9
Bus 4,074,296 44.5 50.8 16.8 25.6 1,677.2 12.5
Custer 4,018,296 39.0 41.3 16.6 20.8 1,849.9 10.3
Auto 2,301,759 22.3 22.3 9.5 11.2 2,997.3 9.7
Combi 2,272,193 18.6 26.0 9.4 13.1 1,032.6 11.5
Mototaxi 2,114,407 4.2 4.2 8.7 2.1 198.5 2.0
Tren 775,325 10.6 14.6 3.2 7.3 0.6 18.8
Taxi 623,171 4.4 4.4 2.6 2.2 813.8 7.1
Colectivo 546,881 5.3 5.3 2.3 2.7 359.8 9.8
Corredor 521,972 6.5 6.9 2.2 3.5 244.1 13.2
Taxi_por_aplicacion 415,550 3.5 3.5 1.7 1.8 639.1 8.4
Metropolitano 329,369 4.7 6.5 1.4 3.3 71.0 19.8
Moto 323,096 3.4 3.4 1.3 1.7 269.2 10.6
Otro 130,953 1.0 1.0 0.5 0.5 95.9 7.3
Alimentador 119,603 2.4 2.7 0.5 1.4 91.9 22.6
Bicicleta 106,556 0.4 0.4 0.4 0.2 0.0 4.2
Patineta 43,669 0.3 0.3 0.2 0.1 0.0 6.4
Camion pequeno 4,907 0.1 0.1 0.0 0.0 24.6 12.3
Camion 3,291 0.1 0.1 0.0 0.0 31.6 23.6
Trailer 1,757 0.1 0.1 0.0 0.0 51.1 42.9
TOTAL 24,227,924 176.2 198.7 100.0 99.9 10,448.2 8.2 
#write.xlsx(estadisticas, "resumen_distancias_emisiones_matriz_con_caminata_previa.xlsx")
estadisticas <- viajes_dist %>% 
  group_by(modo_principal_comparado) %>% 
  summarize(numero_viajes = round(sum(fexp_viaje_nuevo),0), 
            mpkt = round(sum(fexp_viaje_nuevo*Distancia/1000000000),1),
            mpkt_incl_caminata = round(sum(fexp_viaje_nuevo*Distancia_incl_caminata/1000000000),1),
            reparto_modal = round(100*sum(fexp_viaje_nuevo)/sum(viajes_dist$fexp_viaje_nuevo),1), 
            reparto_modal_pkt = round(100*sum(fexp_viaje_nuevo*Distancia_incl_caminata)/sum(viajes_dist$fexp_viaje_nuevo*viajes_dist$Distancia_incl_caminata),1), 
            `CO2-eq`= round(sum(fexp_viaje_nuevo*`CO2-eq`/1000000),1),
            dist_por_viaj = round(sum(fexp_viaje_nuevo*Distancia_incl_caminata/1000)/sum(fexp_viaje_nuevo),1))
  
estadisticas <- estadisticas[order(estadisticas$numero_viajes, decreasing = TRUE),]

total <- data.frame("TOTAL", 
                   sum(estadisticas$numero_viajes),
                   sum(estadisticas$mpkt),
                   sum(estadisticas$mpkt_incl_caminata),
                   sum(estadisticas$reparto_modal),
                   sum(estadisticas$reparto_modal_pkt),
                   sum(estadisticas$`CO2-eq`),
                   round(sum(estadisticas$numero_viajes*estadisticas$dist_por_viaj)/sum(estadisticas$numero_viajes),1))

colnames(total) <- colnames(estadisticas)

estadisticas <- rbind(estadisticas, total)

colnames(estadisticas) <- c("Modo principal","Número de viajes","Millón PKT (solo tramo principal)", "Millón PKT (total incluyendo caminata inicial y final)","% Número de viajes","% PKT","Emisiones GEI (tCO2-eq)", "Distancia media por viaje (km)")

kable(estadisticas, caption = "Estadísticas básicas", format.args = list(big.mark = ","))
Estadísticas básicas
Modo principal Número de viajes Millón PKT (solo tramo principal) Millón PKT (total incluyendo caminata inicial y final) % Número de viajes % PKT Emisiones GEI (tCO2-eq) Distancia media por viaje (km)
Transporte público 9,838,862 107.7 122.8 40.6 61.8 3,934.7 12.5
A pie 5,500,873 4.8 4.8 22.7 2.4 0.0 0.9
Transporte informal 4,933,480 28.1 35.6 20.4 17.9 1,590.9 7.2
Auto 2,301,759 22.3 22.3 9.5 11.2 2,997.3 9.7
Taxi / Carro por aplicación 1,038,721 7.9 7.9 4.3 4.0 1,452.9 7.6
Moto 323,096 3.4 3.4 1.3 1.7 269.2 10.6
Otro 184,578 1.4 1.4 0.8 0.7 203.2 7.8
Bicicleta 106,556 0.4 0.4 0.4 0.2 0.0 4.2
TOTAL 24,227,925 176.0 198.6 100.0 99.9 10,448.2 8.2 
#write.xlsx(estadisticas, "comparado_distancias_emisiones_matriz_con_caminata_previa.xlsx")
estadisticas <- viajes_dist %>% 
  group_by(modo_principal_LATS) %>% 
  summarize(numero_viajes = round(sum(fexp_viaje_nuevo),0), 
            mpkt = round(sum(fexp_viaje_nuevo*Distancia/1000000000),1),
            mpkt_incl_caminata = round(sum(fexp_viaje_nuevo*Distancia_incl_caminata/1000000000),1),
            reparto_modal = round(100*sum(fexp_viaje_nuevo)/sum(viajes_dist$fexp_viaje_nuevo),1), 
            reparto_modal_pkt = round(100*sum(fexp_viaje_nuevo*Distancia_incl_caminata)/sum(viajes_dist$fexp_viaje_nuevo*viajes_dist$Distancia_incl_caminata),1), 
            `CO2-eq`= round(sum(fexp_viaje_nuevo*`CO2-eq`/1000000),1),
            dist_por_viaj = round(sum(fexp_viaje_nuevo*Distancia_incl_caminata/1000)/sum(fexp_viaje_nuevo),1))
  
estadisticas <- estadisticas[order(estadisticas$numero_viajes, decreasing = TRUE),]

total <- data.frame("TOTAL", 
                   sum(estadisticas$numero_viajes),
                   sum(estadisticas$mpkt),
                   sum(estadisticas$mpkt_incl_caminata),
                   sum(estadisticas$reparto_modal),
                   sum(estadisticas$reparto_modal_pkt),
                   sum(estadisticas$`CO2-eq`),
                   round(sum(estadisticas$numero_viajes*estadisticas$dist_por_viaj)/sum(estadisticas$numero_viajes),1))

colnames(total) <- colnames(estadisticas)

estadisticas <- rbind(estadisticas, total)

colnames(estadisticas) <- c("Modo principal","Número de viajes","Millón PKT (solo tramo principal)", "Millón PKT (total incluyendo caminata inicial y final)","% Número de viajes","% PKT","Emisiones GEI (tCO2-eq)", "Distancia media por viaje (km)")

kable(estadisticas, caption = "Estadísticas básicas", format.args = list(big.mark = ","))
Estadísticas básicas
Modo principal Número de viajes Millón PKT (solo tramo principal) Millón PKT (total incluyendo caminata inicial y final) % Número de viajes % PKT Emisiones GEI (tCO2-eq) Distancia media por viaje (km)
Regular bus 11,006,360 110.9 127.8 45.4 64.3 4,895.8 11.6
Walking 5,500,873 4.8 4.8 22.7 2.4 0.0 0.9
Paratransit 2,661,288 9.5 9.5 11.0 4.8 558.3 3.6
Private car 2,301,759 22.3 22.3 9.5 11.2 2,997.3 9.7
Taxi 1,038,721 7.9 7.9 4.3 4.0 1,452.9 7.6
Rapid transit 775,325 10.6 14.6 3.2 7.3 0.6 18.8
BRT 329,369 4.7 6.5 1.4 3.3 71.0 19.8
Moto 323,096 3.4 3.4 1.3 1.7 269.2 10.6
Other 184,578 1.4 1.4 0.8 0.7 203.2 7.8
Bike 106,556 0.4 0.4 0.4 0.2 0.0 4.2
TOTAL 24,227,925 175.9 198.6 100.0 99.9 10,448.3 8.2 
#write.xlsx(estadisticas, "LATS_distancias_emisiones_matriz_con_caminata_previa.xlsx")

Cartografía en Español

Finalmente, producimos mapas por distrito.

limites_lima <- st_make_valid(limites_lima)
centroides_limites_lima <- st_centroid(limites_lima) %>%
  mutate(lon = st_coordinates(.)[,1],
         lat = st_coordinates(.)[,2])

conos_lima <- st_make_valid(conos_lima)
centroides_conos_lima <- st_centroid(conos_lima) %>%
  mutate(lon = c(-77.12,-77.09,-76.75,-77.01,-76.81),
         lat = c(-12.06,-12.13,-12.03,-11.80,-12.17))

centroides_limites_lima$DISPLAY <- centroides_limites_lima$DISTRNOMBR
centroides_limites_lima$DISPLAY[!(centroides_limites_lima$DISTRNOMBR %in% c("VENTANILLA", "SAN JUAN DE LURIGANCHO", "LURIGANCHO (CHOSICA)", "LURIN", "VILLA EL SALVADOR", "LA MOLINA", "SAN BORJA", "MIRAFLORES", "PACHACAMAC", "SANTIAGO DE SURCO", "SANTA ANITA", "LA VICTORIA", "MAGDALENA DEL MAR", "LA PUNTA", "ATE", "INDEPENDENCIA", "CALLAO", "COMAS", "LIMA", "LA PERLA", "CHORRILLOS", "LA PUNTA", "SAN ISIDRO"))] <- ""
centroides_limites_lima$DISPLAY[centroides_limites_lima$DISTRNOMBR == "SAN JUAN DE LURIGANCHO"] <-"SAN JUAN DE \nLURIGANCHO"
centroides_limites_lima$DISPLAY[centroides_limites_lima$DISTRNOMBR == "SAN BORJA"] <-"SAN \nBORJA"
centroides_limites_lima$DISPLAY[centroides_limites_lima$DISTRNOMBR == "SAN ISIDRO"] <-"SAN \nISIDRO"
centroides_limites_lima$DISPLAY[centroides_limites_lima$DISTRNOMBR == "MAGDALENA DEL MAR"] <-"MAGDALENA \nDEL MAR"
centroides_limites_lima$lon[centroides_limites_lima$DISTRNOMBR == "MAGDALENA DEL MAR"] <- -77.09
centroides_limites_lima$lon[centroides_limites_lima$DISTRNOMBR == "SAN BORJA"] <- -76.97
centroides_limites_lima$lon[centroides_limites_lima$DISTRNOMBR == "ATE"] <- -76.84
centroides_limites_lima$lon[centroides_limites_lima$DISTRNOMBR == "PACHACAMAC"] <- -76.87
centroides_limites_lima$lon[centroides_limites_lima$DISTRNOMBR == "VENTANILLA"] <- -77.12
centroides_limites_lima$lon[centroides_limites_lima$DISTRNOMBR == "SAN ISIDRO"] <- -77.025
centroides_limites_lima$lon[centroides_limites_lima$DISTRNOMBR == "LA PUNTA"] <- -77.18
centroides_limites_lima$lat[centroides_limites_lima$DISTRNOMBR == "LA PUNTA"] <- -12.075
centroides_limites_lima$lat[centroides_limites_lima$DISTRNOMBR == "VENTANILLA"] <- -11.91
centroides_limites_lima$lat[centroides_limites_lima$DISTRNOMBR == "LA MOLINA"] <- -12.11
centroides_limites_lima$lat[centroides_limites_lima$DISTRNOMBR == "SANTIAGO DE SURCO"] <- -12.15
centroides_conos_lima$lat[centroides_conos_lima$cono == "Callao"] <- -12.04
Dist_distrito <- viajes_dist[,c("distrito", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_viaje_nuevo")]%>%
  left_join(Distritos[,c("distrito", "pop_est_2024")], by = "distrito") %>%
  group_by(distrito) %>%
  summarize(dist = sum(fexp_viaje_nuevo*Distancia_incl_caminata/1000),
            pop = unique(pop_est_2024),
            dist_capita = dist/pop)

Distritos_no_na <- Dist_distrito[,"distrito"]

Distritos <- Distritos %>% left_join(Dist_distrito[,c("distrito", "dist", "dist_capita")], by = "distrito")
#Colores
bks_total <- round(as.numeric(quantile(Distritos$dist_capita, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),3)
bks <- c(bks_total,max(Distritos$dist_capita))

Distritos$binned_dist_capita <- cut(Distritos$dist_capita,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = Distritos, color = NA, aes(fill = binned_dist_capita))+
  geom_sf(data = limites_lima, fill = NA, aes(color = "Límite de distrito"))+
  geom_text(data= centroides_limites_lima,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 3.5, check_overlap = FALSE) +
  geom_sf(data = conos_lima, fill = NA, aes(color = "Área interdistrital"))+
  geom_text(data= centroides_conos_lima,aes(x=lon, y=lat, label=cono),
            color = "black", fontface = "bold", size = 6, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(6,"Purples")) +
  scale_color_manual(values = c("black", "grey"))+   
  coord_sf(xlim = lon_bounds_Lima, ylim = lat_bounds_Lima, datum = NA)+
  labs(title = "", color = "", fill = "Distancia cotidiana promedio per \ncápita todos modos incluidos (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, ATU, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 91),
        legend.text.align = 1)+
  guides(color = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "bl", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Dist_distrito <- viajes_dist[viajes_dist$modo_principal_comparado == "A pie",c("distrito", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_viaje_nuevo")]%>%
  left_join(Distritos[,c("distrito", "pop_est_2024")], by = "distrito") %>%
  group_by(distrito) %>%
  summarize(dist_walk = sum(fexp_viaje_nuevo*Distancia_incl_caminata/1000),
            pop = unique(pop_est_2024),
            walk_intensity = dist_walk/pop)

Distritos <- Distritos %>% left_join(Dist_distrito[,c("distrito", "dist_walk", "walk_intensity")], by = "distrito")
#Colores
bks_walk <- round(as.numeric(quantile(Distritos$walk_intensity, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),10)
bks <- c(bks_walk,max(Distritos$walk_intensity))

Distritos$binned_walk_intensity <- cut(Distritos$walk_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = Distritos, color = NA, aes(fill = binned_walk_intensity))+
  geom_sf(data = limites_lima, fill = NA, aes(color = "Límite de distrito"))+
  #geom_text(data= centroides_limites_lima,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 3.5, check_overlap = FALSE) +
  geom_sf(data = conos_lima, fill = NA, aes(color = "Área interdistrital"))+
  #geom_text(data= centroides_conos_lima,aes(x=lon, y=lat, label=cono), color = "black", fontface = "bold", size = 6, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(6,"Reds")) +
  scale_color_manual(values = c("black", "grey"))+   
  coord_sf(xlim = lon_bounds_Lima, ylim = lat_bounds_Lima, datum = NA)+
  labs(title = "", color = "", fill = "Distancia cotidiana promedio \nper cápita caminando (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, ATU, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 91),
        legend.text.align = 1)+
  guides(color = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "bl", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Dist_distrito <- viajes_dist[viajes_dist$modo_principal_comparado == "Auto",c("distrito", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_viaje_nuevo")]%>%
  left_join(Distritos[,c("distrito", "pop_est_2024")], by = "distrito") %>%
  group_by(distrito) %>%
  summarize(dist_auto = sum(fexp_viaje_nuevo*Distancia_incl_caminata/1000),
            pop = unique(pop_est_2024),
            auto_intensity = dist_auto/pop)

Distritos <- Distritos %>% left_join(Dist_distrito[,c("distrito", "dist_auto", "auto_intensity")], by = "distrito")
#Colores
bks_car <- round(as.numeric(quantile(Distritos$auto_intensity, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),5)
bks <- c(bks_car,max(Distritos$auto_intensity))

Distritos$binned_auto_intensity <- cut(Distritos$auto_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = Distritos, color = NA, aes(fill = binned_auto_intensity))+
  geom_sf(data = limites_lima, fill = NA, aes(color = "Límite de distrito"))+
  #geom_text(data= centroides_limites_lima,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 3.5, check_overlap = FALSE) +
  geom_sf(data = conos_lima, fill = NA, aes(color = "Área interdistrital"))+
  #geom_text(data= centroides_conos_lima,aes(x=lon, y=lat, label=cono), color = "black", fontface = "bold", size = 6, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(6,"RdPu")) +
  scale_color_manual(values = c("black", "grey"))+   
  coord_sf(xlim = lon_bounds_Lima, ylim = lat_bounds_Lima, datum = NA)+
  labs(title = "", color = "", fill = "Distancia cotidiana promedio per \ncápita en auto particular (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, ATU, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 91),
        legend.text.align = 1)+
  guides(color = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "bl", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Dist_distrito <- viajes_dist[viajes_dist$modo_principal_comparado == "Transporte público",c("distrito", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_viaje_nuevo")]%>%
  left_join(Distritos[,c("distrito", "pop_est_2024")], by = "distrito") %>%
  group_by(distrito) %>%
  summarize(dist_transit = sum(fexp_viaje_nuevo*Distancia_incl_caminata/1000),
            pop = unique(pop_est_2024),
            transit_intensity = dist_transit/pop)

Distritos <- Distritos %>% left_join(Dist_distrito[,c("distrito", "dist_transit", "transit_intensity")], by = "distrito")
#Colores
bks_transit <- round(as.numeric(quantile(Distritos$transit_intensity, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- c(bks_transit,max(Distritos$transit_intensity))

Distritos$binned_transit_intensity <- cut(Distritos$transit_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = Distritos, color = NA, aes(fill = binned_transit_intensity))+
  geom_sf(data = limites_lima, fill = NA, aes(color = "Límite de distrito"))+
  #geom_text(data= centroides_limites_lima,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 3.5, check_overlap = FALSE) +
  geom_sf(data = conos_lima, fill = NA, aes(color = "Área interdistrital"))+
  #geom_text(data= centroides_conos_lima,aes(x=lon, y=lat, label=cono), color = "black", fontface = "bold", size = 6, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(6,"Greens")) +
  scale_color_manual(values = c("black", "grey"))+   
  coord_sf(xlim = lon_bounds_Lima, ylim = lat_bounds_Lima, datum = NA)+
  labs(title = "", color = "", fill = "Distancia cotidiana promedio per \ncápita en transporte público (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, ATU, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 91),
        legend.text.align = 1)+
  guides(color = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "bl", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Dist_distrito <- viajes_dist[,c("distrito", "modo_principal_comparado", "CO2-eq", "fexp_viaje_nuevo")]%>%
  left_join(Distritos[,c("distrito", "pop_est_2024")], by = "distrito") %>%
  group_by(distrito) %>%
  summarize(co2 = sum(fexp_viaje_nuevo*`CO2-eq`/1000000),
            pop = unique(pop_est_2024),
            co2_capita = 1000*co2/pop)

Dist_distrito <- Distritos_no_na %>% left_join(Dist_distrito, by = "distrito")
Dist_distrito$co2_capita[is.na(Dist_distrito$co2_capita)] <-0

Distritos <- Distritos %>% left_join(Dist_distrito[,c("distrito", "co2", "co2_capita")], by = "distrito")
#Colores
bks_co2 <- round(as.numeric(quantile(Distritos$co2_capita, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- c(bks_co2,max(Distritos$co2_capita))

Distritos$binned_co2_capita <- cut(Distritos$co2_capita,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = Distritos, color = NA, aes(fill = binned_co2_capita))+
  geom_sf(data = limites_lima, fill = NA, aes(color = "Límite de distrito"))+
  #geom_text(data= centroides_limites_lima,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 3.5, check_overlap = FALSE) +
  geom_sf(data = conos_lima, fill = NA, aes(color = "Área interdistrital"))+
  #geom_text(data= centroides_conos_lima,aes(x=lon, y=lat, label=cono), color = "black", fontface = "bold", size = 6, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(6,"Oranges")) +
  scale_color_manual(values = c("black", "grey"))+   
  coord_sf(xlim = lon_bounds_Lima, ylim = lat_bounds_Lima, datum = NA)+
  labs(title = "", color = "", fill = "Emisiones cotidianas promedio \nper cápita (kg CO2-eq)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, ATU, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 91),
        legend.text.align = 1)+
  guides(color = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "bl", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Cartography in English

#Colores
bks_total <- round(as.numeric(quantile(Distritos$dist_capita, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),3)
bks <- c(bks_total,max(Distritos$dist_capita))

Distritos$binned_dist_capita <- cut(Distritos$dist_capita,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = Distritos, color = NA, aes(fill = binned_dist_capita))+
  geom_sf(data = limites_lima, fill = NA, aes(color = "District boundary"))+
  geom_text(data= centroides_limites_lima,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 3.5, check_overlap = FALSE) +
  geom_sf(data = conos_lima, fill = NA, aes(color = "Interdistrict area"))+
  geom_text(data= centroides_conos_lima,aes(x=lon, y=lat, label=cono),
            color = "black", fontface = "bold", size = 6, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(5,"Purples")) +
  scale_color_manual(values = c("black", "grey"))+   
  coord_sf(xlim = lon_bounds_Lima, ylim = lat_bounds_Lima, datum = NA)+
  labs(title = "", color = "", fill = "Average daily distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Lima & Callao Urban Transport Authority (ATU), 2023 \nDiscretization in quintiles", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
        legend.text.align = 1)+
  guides(color = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "bl", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

#Colores
bks_walk <- round(as.numeric(quantile(Distritos$walk_intensity, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),10)
bks <- c(bks_walk,max(Distritos$walk_intensity))

Distritos$binned_walk_intensity <- cut(Distritos$walk_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = Distritos, color = NA, aes(fill = binned_walk_intensity))+
  geom_sf(data = limites_lima, fill = NA, aes(color = "District boundary"))+
  geom_sf(data = conos_lima, fill = NA, aes(color = "Interdistrict area"))+
  scale_fill_manual(values = brewer.pal(5,"Reds")) +
  scale_color_manual(values = c("black", "grey"))+   
  coord_sf(xlim = lon_bounds_Lima, ylim = lat_bounds_Lima, datum = NA)+
  labs(title = "", color = "", fill = "Average daily walking distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Lima & Callao Urban Transport Authority (ATU), 2023 \nDiscretization in quintiles", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
        legend.text.align = 1)+
  guides(color = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "bl", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

#Colores
bks_car <- round(as.numeric(quantile(Distritos$auto_intensity, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),5)
bks <- c(bks_car,max(Distritos$auto_intensity))

Distritos$binned_auto_intensity <- cut(Distritos$auto_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = Distritos, color = NA, aes(fill = binned_auto_intensity))+
  geom_sf(data = limites_lima, fill = NA, aes(color = "District boundary"))+
  geom_sf(data = conos_lima, fill = NA, aes(color = "Interdistrict area"))+
  scale_fill_manual(values = brewer.pal(5,"RdPu")) +
  scale_color_manual(values = c("black", "grey"))+   
  coord_sf(xlim = lon_bounds_Lima, ylim = lat_bounds_Lima, datum = NA)+
  labs(title = "", color = "", fill = "Average daily private car distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Lima & Callao Urban Transport Authority (ATU), 2023 \nDiscretization in quintiles", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
        legend.text.align = 1)+
  guides(color = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "bl", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

#Colores
bks_transit <- round(as.numeric(quantile(Distritos$transit_intensity, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),2)
bks <- c(bks_transit,max(Distritos$transit_intensity))

Distritos$binned_transit_intensity <- cut(Distritos$transit_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = Distritos, color = NA, aes(fill = binned_transit_intensity))+
  geom_sf(data = limites_lima, fill = NA, aes(color = "District boundary"))+
  geom_sf(data = conos_lima, fill = NA, aes(color = "Interdistrict area"))+
  scale_fill_manual(values = brewer.pal(5,"Greens")) +
  scale_color_manual(values = c("black", "grey"))+   
  coord_sf(xlim = lon_bounds_Lima, ylim = lat_bounds_Lima, datum = NA)+
  labs(title = "", color = "", fill = "Average daily public transport distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Lima & Callao Urban Transport Authority (ATU), 2023 \nDiscretization in quintiles", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
        legend.text.align = 1)+
  guides(color = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "bl", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

#Colores
bks_co2 <- round(as.numeric(quantile(Distritos$co2_capita, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),2)
bks <- c(bks_co2,max(Distritos$co2_capita))

Distritos$binned_co2_capita <- cut(Distritos$co2_capita,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = Distritos, color = NA, aes(fill = binned_co2_capita))+
  geom_sf(data = limites_lima, fill = NA, aes(color = "District boundary"))+
  geom_sf(data = conos_lima, fill = NA, aes(color = "Interdistrict area"))+
  scale_fill_manual(values = brewer.pal(5,"Oranges")) +
  scale_color_manual(values = c("black", "grey"))+   
  coord_sf(xlim = lon_bounds_Lima, ylim = lat_bounds_Lima, datum = NA)+
  labs(title = "", color = "", fill = "Average daily GHG emissions \nper capita (kg CO2-eq)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Lima & Callao Urban Transport Authority (ATU), 2023 \nDiscretization in quintiles", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
        legend.text.align = 1)+
  guides(color = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "bl", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Correlaciones

Distritos_ICS <- read.xlsx("Distritos_ICS.xlsx")

Distritos <- Distritos %>% left_join(Distritos_ICS[, c("distrit", "propICS12")], by = c("distrito" = "distrit"))

correlaciones <- rbind(data.frame(cor.test(Distritos$propICS12, Distritos$auto_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(Distritos$walk_intensity, Distritos$auto_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(Distritos$walk_intensity, Distritos$transit_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(Distritos$propICS12, Distritos$walk_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(Distritos$propICS12, Distritos$transit_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(Distritos$co2_capita, Distritos$auto_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(Distritos$propICS12, Distritos$co2_capita, method = "pearson")[c(4,3)]))

rownames(correlaciones) <- c("auto_ICS",
                             "walk_auto",
                             "walk_transit",
                             "ICS_walk",
                             "ICS_transit",
                             "GES_auto",
                             "ICS_GES")

kable(correlaciones, caption = "Correlaciones", format.args = list(big.mark = ","))
Correlaciones
estimate p.value
auto_ICS -0.4003704 0.0039639
walk_auto -0.1985772 0.1668337
walk_transit -0.0598660 0.6796209
ICS_walk 0.1332195 0.3563731
ICS_transit 0.5807414 0.0000098
GES_auto 0.7163892 0.0000000
ICS_GES -0.0112065 0.9384324 
correlaciones <- cbind(rownames(correlaciones), correlaciones)

write.xlsx(correlaciones, "correlaciones.xlsx", overwrite = TRUE)

Palma ratio

Distritos_top_10 <- Distritos[order(Distritos$propICS12, decreasing = FALSE),]
select <-as.integer(0.1*nrow(Distritos_top_10))
Distritos_top_10 <- Distritos_top_10[1:select,]

Distritos_bottom_40 <- Distritos[order(Distritos$propICS12, decreasing = TRUE),]
select <-as.integer(0.4*nrow(Distritos_bottom_40))
Distritos_bottom_40 <- Distritos_bottom_40[1:select,]

palma <- rbind(data.frame(top_10 = weighted.mean(x = Distritos_top_10$dist_capita, w = Distritos_top_10$pop_est_2024),
                   bottom_40 = weighted.mean(x = Distritos_bottom_40$dist_capita, w = Distritos_bottom_40$pop_est_2024)),
               data.frame(top_10 = weighted.mean(x = Distritos_top_10$auto_intensity, w = Distritos_top_10$pop_est_2024),
                   bottom_40 = weighted.mean(x = Distritos_bottom_40$auto_intensity, w = Distritos_bottom_40$pop_est_2024)),
               data.frame(top_10 = weighted.mean(x = Distritos_top_10$walk_intensity, w = Distritos_top_10$pop_est_2024),
                   bottom_40 = weighted.mean(x = Distritos_bottom_40$walk_intensity, w = Distritos_bottom_40$pop_est_2024)),
               data.frame(top_10 = weighted.mean(x = Distritos_top_10$transit_intensity, w = Distritos_top_10$pop_est_2024),
                   bottom_40 = weighted.mean(x = Distritos_bottom_40$transit_intensity, w = Distritos_bottom_40$pop_est_2024)),
               data.frame(top_10 = weighted.mean(x = Distritos_top_10$co2_capita, w = Distritos_top_10$pop_est_2024),
                   bottom_40 = weighted.mean(x = Distritos_bottom_40$co2_capita, w = Distritos_bottom_40$pop_est_2024)))


palma$palma <- palma$top_10/palma$bottom_40

rownames(palma) <- c("dist_capita",
                     "auto_intensity",
                     "walk_intensity",
                     "transit_intensity",
                     "co2_capita")

kable(palma, caption = "Palma ratios", format.args = list(big.mark = ","))
Palma ratios
top_10 bottom_40 palma
dist_capita 13.8238132 19.7752379 0.6990466
auto_intensity 5.9482566 1.2994395 4.5775557
walk_intensity 0.3645133 0.4139389 0.8805968
transit_intensity 5.0871109 12.7297614 0.3996234
co2_capita 1.2667923 0.8709589 1.4544801 
palma <- cbind(rownames(palma), palma)

write.xlsx(palma, "palma.xlsx", overwrite = TRUE)

Valores promediados de los indicadores

promedios <- data.frame(dist_capita = sum(Distritos$dist)/sum(Distritos$pop_est_2024),
                        walk_intensity = sum(Distritos$dist_walk)/sum(Distritos$pop_est_2024),
                        auto_intensity = sum(Distritos$dist_auto)/sum(Distritos$pop_est_2024),
                        transit_intensity = sum(Distritos$dist_transit)/sum(Distritos$pop_est_2024),
                        ges_capita = 1000*sum(Distritos$co2)/sum(Distritos$pop_est_2024))

kable(promedios, caption = "Valores promediados de los indicadores", format.args = list(big.mark = ","))
Valores promediados de los indicadores
dist_capita walk_intensity auto_intensity transit_intensity ges_capita
17.67728 0.4294056 1.981313 10.92735 0.9294103